function [p,w,ZMat,ZstarMat, LMat, YMat, NMat, IMat, VMat, ZbarMat, Cratio] = Counterfactuals(pm,bc,bp)

pm.Xmax = 10;
gpts = pm.Xmax./0.0001;

options = optimoptions('fsolve','Display','off','Algorithm','trust-region','FunctionTolerance',1e-8,'StepTolerance',1e-8);
p =pm.price;
w=pm.wage;
tol = 10; iter = 1; 
%conv = 0;
%ze0 =logninv(1-(Data.N_f+Data.N_i)./pm.M,0,pm.sigma);
%zr0 = logninv(1-Data.N_f./pm.M,0,pm.sigma);

pm.bc = bc;
pm.bp = bp;
Ee_i= zeros(pm.I,1);
for r = 1:pm.I
    Ee_i(r)= quadgk(@(e) e.^(pm.theta(r)./(pm.theta(r)-pm.rho)).*lognpdf(e,0,pm.se),-Inf,Inf);
end

%if informal == 1 
%    fprintf('\n \n There is informal sector \n \n');
    while tol>pm.mintolp 
            %wf = pm.B.*w;

            % ze= fsolve(@(x) pi_i(pm,x,p,w,0) -p.*pm.E_i,ze0,options);
            % zr = fsolve(@(x) pi_f(pm,x,p,wf,pm.t).*Ee - pi_i(pm,x,p,w,0)-p.*pm.E_r,zr0,options);
                        
            % CEF = fsolve(@(x) pi_f(pm,x,p,wf,pm.t).*Ee -p.*(pm.E_i+pm.E_r),ones(pm.I,1),options);
            % zr = max(CEF,WEF);

            % 
             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             % CHOICES FOR FIRMS
             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %     gpts = 10
                X = repmat(linspace(0,pm.Xmax,gpts),pm.I,1);

                PiMat = zeros(pm.I,gpts,3);
                PiMat(:,:,1)= 0;
                
                for r = 1:pm.I
                    PiMat(r,:,2) = integral(@(e) (pi_i(pm,X(r,:),e,p(r),w,r) - p(r).*pm.E_i(r)) ...
                         .*lognpdf(e,0,pm.se),-Inf,Inf,'ArrayValued',true);
                    PiMat(r,:,3) = integral(@(e) (pi_f(pm,X(r,:),e,p(r),w,pm.t(r),r) - p(r).*(pm.E_i(r)+pm.E_r(r))) ...
                         .*lognpdf(e,0,pm.se),-Inf,Inf,'ArrayValued',true);
                end

                Pimax = max(PiMat,[],3);
                ChoiceMat = PiMat == Pimax;

                Xe= X.*ChoiceMat(:,:,2);
                Xe(Xe==0)=NaN;
                ze = min(Xe,[],2);

                Xf= X.*ChoiceMat(:,:,3);
                Xf(Xf==0)=NaN;
                zr = min(Xf,[],2);
                ze(isnan(ze)) = zr(isnan(ze));
                zr(isnan(zr)) = ze(isnan(zr));
            %}

            seedha = ze<zr;
            %Total profits in the formal and informal sectors 
            L_i = zeros(pm.I,1);
            L_f = zeros(pm.I,1);
            L_c = zeros(pm.I,1);
            L_p = zeros(pm.I,1);
            wfL_f = zeros(pm.I,1);
            Y_i = zeros(pm.I,1);
            Y_f = zeros(pm.I,1);
            Pi_i = zeros(pm.I,1);
            Pi_f = zeros(pm.I,1);
        
      
        for r = 1:pm.I
            L_i(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x)  l_i(pm,x,e,p(r),w,r)...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf);
            
            L_f(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x)  l_f(pm,x,e,p(r),w,pm.t(r),r)...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);

            wfL_f(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x)  ...
                    wf(pm,x,e,w,r).*l_f(pm,x,e,p(r),w,pm.t(r),r)...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);

            Y_i(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x) x.*e.* l_i(pm,x,e,p(r),w,r).^pm.rho...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf);
            
            Y_f(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x) x.*e.*  l_f(pm,x,e,p(r),w,pm.t(r),r).^pm.rho...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);

            Pi_i(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x) (pi_i(pm,x,e,p(r),w,r) - p(r).*pm.E_i(r)) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf);
            
            Pi_f(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x) (pi_f(pm,x,e,p(r),w,pm.t(r),r) - p(r).*(pm.E_i(r)+pm.E_r(r)))...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);

            L_c(r) = integral(@(e) integral(@(x) ...
                    (w.*pm.bc(r)./((x.*e).^pm.phi.*wf(pm,x,e,w,r))).^(-pm.nu) .*L_f(r) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),Inf,'ArrayValued',true),-Inf,Inf);

            L_p(r) = integral(@(e) integral(@(x) ...
                    (w.*pm.bp(r)./wf(pm,x,e,w,r)).^(-pm.nu) .*L_f(r) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);
        end

        
            %Fixed costs and Taxes paid
            FC = p.*pm.E_i.*pm.M.*(1-logncdf(min(ze,zr),0,pm.sigma)) + p.*pm.E_r.*pm.M.*(1-logncdf(zr,0,pm.sigma));
            Tpaid = pm.t.*p.*Y_f;

%            income = w.*sum(L_i + L_c + L_p,'omitnan') + sum(Pi_i + Pi_f + Tpaid,'omitnan');
        income = sum(w.*L_i + wfL_f,'omitnan') + sum(Pi_i + Pi_f + Tpaid,'omitnan');
%            income = w.*sum(L_i + bc.*L_c + bp.*L_p,'omitnan') + sum(Pi_i + Pi_f + Tpaid,'omitnan');
            pnew = (pm.kappa.*income + FC)./ (Y_i+Y_f);
            pnew(pnew>1000)=1000;
            Pnew = exp(sum(pm.kappa.*log(pnew./pm.kappa)));
            pnew = pnew./Pnew;

            frac = sum(L_i+L_c+L_p,'omitnan')./pm.L;
            tolp = abs((p-pnew))./p;
            tolw = abs(frac-1);
            a= 0.8;
            p = a.*p +(1-a).*pnew;
            w = a.*w +(1-a).*frac.*w;
    %        p(1) = pm.alpha(1).*exp(-1./pm.alpha(1) .*sum(pm.alpha(2:pm.I).*log(p(2:pm.I)./pm.alpha(2:pm.I))));
            tol = max(max(tolp),tolw);
            conv = sum(tolp<=pm.mintolp)+sum(tolw<=pm.mintolp) - 1;
            fprintf('Iter = %d; Converged = %d out of %d; tol = %0.4f \n',iter,conv,pm.I,tol);
    %        [tolw;tolp]
            iter = iter+1;
    end


ZstarMat = [ze zr];
NMat= pm.M.*[logncdf(zr,0,pm.sigma)-logncdf(ze,0,pm.sigma) ...
                1-logncdf(zr,0,pm.sigma)];
Cratio = zeros(pm.I,1);

ZMat = zeros(pm.I,2);
ZbarMat = zeros(pm.I,1);
for r = 1:pm.I
    ZbarMat(r,:) = pm.M(r).*integral(@(x) x.*lognpdf(x,0,pm.sigma),ze(r),Inf)./(1-logncdf(ze(r),0,pm.sigma));
    ZMat(r,:) = pm.M(r).*[integral(@(x) x.*lognpdf(x,0,pm.sigma),ze(r),zr(r)) ...
                integral(@(x) x.*lognpdf(x,0,pm.sigma),zr(r),Inf)];
%    ZMat(r,:) = pm.M(r).*[(logncdf(zr(r),0,pm.sigma)-logncdf(ze(r),0,pm.sigma)) ...
%               (1-logncdf(zr(r),0,pm.sigma))];
    xbar = integral(@(x) x.^pm.phi.*lognpdf(x,0,pm.sigma),zr(r),Inf) ...
                ./(1-logncdf(zr(r),0,pm.sigma));
    Ee= integral(@(e) e.^pm.phi.*lognpdf(e,0,pm.se),-Inf,Inf);
    Cratio(r) = (pm.bc./(xbar.*Ee.*pm.bp)).^(-pm.nu);
end

varz_i = zeros(pm.I,1);
varz_f = zeros(pm.I,1);
for r = 1:pm.I
            lbar_i = integral(@(e) integral(@(x)  log(l_i(pm,x,e,p(r),w,r)) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf) ...
                ./(logncdf(zr(r),0,pm.sigma)-logncdf(ze(r),0,pm.sigma));

            lbar_f = integral(@(e) integral(@(x)  log(l_f(pm,x,e,p(r),w,pm.t(r),r)) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf) ...
                ./(1 - logncdf(zr(r),0,pm.sigma));

            %variance of x
            varz_i(r) = integral(@(e) integral(@(x)  (log(l_i(pm,x,e,p(r),w,r)) - lbar_i).^2 ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf) ...
                ./(logncdf(zr(r),0,pm.sigma)-logncdf(ze(r),0,pm.sigma));

            varz_f(r) = integral(@(e) integral(@(x)  (log(l_f(pm,x,e,p(r),w,pm.t(r),r)) - lbar_f).^2 ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf) ...
                ./(1 - logncdf(zr(r),0,pm.sigma));
end
VMat = [varz_i; varz_f]; 


LMat = [L_i L_c+L_p];
YMat = [Y_i Y_f];
IMat =  w.*sum(L_i+ L_c + L_p,'omitnan') + sum(Pi_i + Pi_f + Tpaid,'omitnan');
%TFP = (sum(Y_i + Y_f))./sum((L_i + L_f).^(pm.rho));
